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On  Multichannel  (Multivariate)  Maximum  Entropy  Spectral  Analysis 
C.H.  Chen  and  Chihsung  Yen 

1.  Introduction 

The  univariate  maximum  entropy  spectral  analysis  has  now  been  well  developed 
and  applied  to  many  defense  research  areas  such  as  radar,  sonar  and  geophysics.  There 
has  been  some  work  done  to  extend  the  maximum  entropy  spectral  analysis  to  multivariate 
case.  Whittle  [1]  and  Robinson  [2] [3]  generalized  the  Levinson-Durbin  recursion  to 
the  multivariate  case  by  fitting  both  foward  and  backward  autoregressions  in  a 
stepwise  fashion.  In  this  thesis,  Burg  [4]  has  mentioned  about  the  multichannel  case. 
However  the  computer  programs  for  both  multichannel  and  multivariate  maximum  entropy 
spectral  analysis  were  only  recently  developed  successfully.  Morf  et.al  [5]  developed 
an  algorithm  for  direct  estimation  of  the  normalized  reflection  coefficients  from  the 
observed  data  for  maximum  entropy  spectral  analysis.  They  also  compared  the  spectral 
estimation  with  the  methods  of  Jones  [6],  Nuttall  [7]  and  Strand  [8],  which  are  more 
of  a  direct  extension  of  Burg's  work  to  the  multichannel  (multivariate)  case.  Burg's 
algorithm  does  not  generalize  directly  since  the  forward  and  backward  autoregression 
matrices  are  not  the  same  in  the  multivariate  case,  and  the  forward  and  backward  one- 
step  prediction  error  covariance  matrices  are  different,  although  they  have  the  same 
determinant.  In  this  report,  the  programs  developed  by  Strand  and  Jones  are  applied 
to  real  multichannel  data  and  imagery  data  in  addition  to  a  set  of  test  signals. 

The  merits  of  these  methods  are  closely  examined.  In  spite  of  programming  complexity 
the  multichannel  and  multivariate  maximum  entropy  spectral  analysis  will  have  increased 
application  as  the  real  data  are  almost  always  gathered  in  several  channels.  Data 
from  several  channels  form  a  vector  for  multivariate  study. 

2.  Brief  Mathematical  Analysis 

Let  x^.x^, ...,xn  denote  n  zero  mean  vectors  of  dimension  d  each.  The  sample 
estimate  of  covariance  sequence  for  lag  j  is 
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where  the  prime  denotes  the  transposed  vector.  The  forward  and  backward  predicting 


autoregressions  of  order  p  are  given,  respectively,  as 
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where  AR*^  and  BRP^  are  d  x  d  matrices,  and  can  be  determined  recursively  [6]  by 


making  use  of  the  estimated  covariance  matrix  in  Eq.  (1) .  The  recursion  starts  with 

s(f)=  s(t)=  R  (3) 

o  o  o 

The  one-step  forward  and  backward  prediction  error  covariance  matrices  are 
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The  forward  and  backward  residuals  are,  respectively ,  / 
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The  recursive  equations  are  then  given  by 

t  =  p+1,, 


e(p)  =  e(P-1>  _  a(p)b(p_1) 
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The  least  squares  estimates  for  the  forward  and  backward  autoregression  matrices  are 

(7) 
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where  U  is  the  sum  of  cross  products  of  forward  and  backward  residuals  at  lag  p, 
nlP 

(8) 
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and  V  and  W  are  estimates  of  (n  -  p)S^  ,  (n  -  p)S^^  respectively, 


v  =  IP  e*p_1)  (e^p_1))' 


(io) 

Although  the  forward  and  backward  autoregression  matrices  and  the  prediction  error 
covariance  matrices  are  different,  the  multivariate  spectra  should  be  identical 
when  calculated  from  the  forward  and  backward  fits  by 


S (f)  =  h[A(f)]_1  S^f)  [A*(f)]_1 


or  by 


where 


S(f)  =  h[B(f) ]_1  S<b)  [B*(f)]_1 
A(f)  =  I  -  l  a[p)  e27tikhf 
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h  is  the  sampling  period  and  *  denotes  complex  conjugate  transpose. 

The  above  approach  based  on  the  work  of  Jones  [6]  does  not  guarantee  stability 
and  does  not  generally  produce  a  non-negative  definite  spectrum  as  has  been  pointed 

J 

out  by  Nuttall  [7],  Subsequently,  Nuttall  [7]  and  Strand  [8]  applied  a  weighted 
arithmetic  mean  error  criterion  in  order  to  provide  model  stability  and  to  ensure 
positive  definite  stationary  spectra.  Another  procedure  suggested  by  Morf,  et.al. 
[5]  that  also  meets  the  spectral  requirements  is  to  compute  the  spectrum  from  the 
normalized  reflection  coefficient  matrix  p.  To  obtain  this  matrix,  W  and  V  are 
factored  by  using  Cholesky  decomposition  into  the  product  of  lower  triangular 
matrices  times  their  transposes.  A  new  recursive  procedure  for  S ^  and  S^b^  by 
using  p  in  place  of  Eqs.  (4)  can  then  be  obtained.  Other  recursive  algorithm  has 
been  proposed  [9]  for  the  solution  of  the  normal  equations  for  both  single  and 
multichannel  data. 


i 


3.  Spectral  Analysis  of  Multichannel  Sinusoids 

To  test  his  multichannel  maximum  entropy  spectral  analysis  program.  Strand  [8] 
defined  a  two-channel  complex  sinusoid  with  frequency  0.0625  Hz.  Each  channel  con¬ 
tains  128  complex  numbers  plus  a  small  amount  of  additive  Gaussian  noise.  Morf,  et. 
al«  employed  the  same  test  signals  and  demonstrated  that  their  multichannel  algorithm 
performs  better  than  that  of  Strand.  The  complete  data  set  is  tabulated  in  Appendix 
I,  In  our  study  the  data  set  is  considered  as  four-channel  real  sinusoids.  Fig.  1 
is  a  plot  of  the  first  two  channels.  The  computer  program  developed  in  our  study 
is  based  on  the  mathematical  analysis  of  the  previous  section  and  a  subroutine  due 
to  Strand.  At  9  lags,  our  spectral  peak  value  is  2100  as  compared  to  2000  at  the 
frequency  of  0.0625  Hz,  obtained  by  Morf.  For  9  lags.  Fig.  2a  is  the  spectral  plot 
in  linear  scale  while  Fig.  2b  uses  the  logarithmic  scale.  In  Fig.  2c  the  upper 
photo  shows  the  phase  difference  while  the  lower  photo  is  the  coherence  function 
between  the  two  channels.  If  all  four  channels  are  considered  simultaneously,  the 
spectral  peaks  will  be  proportionally  reduced  while  the  shapes  of  the  spectra  remain 
almost  unchanged. 

^•comparison  is  also  made  with  the  Jones'  method  by  using  the  computer  program 
listed  in  the  Appendix  of  Ref.  6.  The  optimum  order  (i.e.  the  number  of  lags)  is 
found  to  be  7  for  the  four-channel  test  data  studied.  Fig.  3a  is  the  spectral  plot 
of  channels  1  &  2  in  linear  scale  (left)  and  in  logarithmic  scale  (right).  Similarly 
Fig.  3b  is  the  spectral  plot  of  channels  3  &  4.  Fig.  3c  is  a  plot  of  phase 
differences  between  two-channels.  Fig.  3d  is  a  plot  of  coherence  functions  between 
two-channels.  Although  there  is  little  change  in  the  spectral  shapes,  the  Jones' 
method  provides  much  smaller  spectral  peaks  (Note  that  all  figures  are  normalized 
plots) .  Thus  the  computer  program  given  in  the  Appendix  provides  somewhat  better 
results  than  those  available  in  other  methods. 


Fig.  3a 
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4.  Spectral  Analysis  of  Geophysical  Data 

A  three-channel  geophysical  data  set  is  provided  by  Robinson  [10]  that  includes 
sunspot  numbers,  northern  light  activity  and  earthquake  activity  with  100  data  points 
each.  Fig.  4b  is  the  spectral  plot  of  each  channel  in  linear  scale,  while  Fig.  4c  is 
the  same  plot  in  logarithmic  scale.  Fig.  4d  shows  the  phase  differences  between  two- 
channels.  The  number  of  lags  is  17  which  is  determined  as  the  optimum  order.  It  is 
interesting  to  note  that  spectral  peak  for  the  sunspot  numbers  is  determined 
accurately. 

5.  Image  Segmentation  Via  Multivariate  Spectral  Analysis 

The  use  of  multivariate  autoregressive  analysis  for  texture  classification  and 
segmentation  was  considered  by  Gambotto  and  Gueguen  [11].  Their  segmentation  result 
however  is  much  to  be  desired.  We  have  employed  the  multivariate  autoregression 
estimation  program  of  Jones  [6]  on  the  infrared  image  as  shown  in  Fig.  5a.  A  small 
area  of  64x64  is  considered  in  the  computer  study.  The  subimage  is  processed  by 
small  squares  of  2x10  pixels  each,  making  use  of  the  two-channel  multivariate 
autoregression  analysis.  Then  the  mean-square  prediction  error  (defined  by  Eq.14  of 
Ref.  11)  is  computed  and  compared  with  a  threshold.  The  idea  is  that  the  prediction 
error  tends  to  be  large  for  pixels  near  the  object  boundary  and  it  is  small  for 
homogeneous  region.  Segmentation  is  performed  by  pixel  classification  using  the 
prediction  error.  The  classified  result  superimposed  on  the  binary  image  of  original 
picture  is  shown  in  Fig,  5b  for  orders  4,  5,  6  and  7.  These  results  are  acceptable 
but  not  as  good  as  compared  with  the  segmentation  using  Figher’s  linear  discrimunant 
112],  Furthermore  the  amount  of  computation  required  for  the  multivariate  auto¬ 
regression  analysis  is  also  much  larger.  Further  study  of  image  analysis  using 
multivariate  autoregression  is  much  needed. 


Fig.  4a 
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Appendix  X  Four-channel 
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Appendix  II  Program  MULl  (Multi-channel  Maximum  Entropy  Spectral  Analysis) 
MUL ) 

REAL  ET  <  4 ,  1 28 ) .  BY  ( 4 ,  1 28 ) ,  UU  ( 4 ) 

COMMON  / BLKO/LT,  BT 
COMMON  INDEX 

DEFINE  FILE  7(  128,  8,  U,  LINE ) 

READ <6, 20 )N, NR.  NSPE. NPRS.  NFREQ 
20  FORMAT (516) 

NOUT  =  1 28 
NDAT=128 
DT*1. 

DO  2  K=  1 ,  NDA'T 
LINE  =K 

READ  (7  -MNFDUU 
DO  3  1-1, NP 
ET ( I . K ) -UU ( I ) 

3  BT ( I .  K ) =UU ( I ) 

2  CONTINUE 

END  FILE  7 

DEF I NE  FILE  1(16,  256 ,  U ,  I NDFI X  ) 

CALL  STRAND (NOUT.  NP,  NUAT.  N,  BT,  NFREQ,  NSPE,  NPRS) 

CALL  EXIT 
END 

SUBROUTINE  CPRINT(NP, ND AT ,  19,  P, PP,  CN, CNP,  CF, CB) 

REAL  P  ( 4,  4 ) ,  PP  ( 4,  4  > ,  CN  (4,4),  CNP  (4,4),  CF  ( 4,  4,  10) .  CB  (4,4,  10 ) 

REAL  ET ( 4 , 128), BT ( 4 , 128) 

COMMON  /BI..K0/ET,  BT 
I 9P 1-19+1 
NDATM9=NDAT~ I 9 
PRINT  2,  19,  I9P1 

2  FORM  AT  ( 1  HI, ///•'  RESULTS  BELOW  FOR'-,  15,  •'  LAGS,  I.  E.  MATRIX 
1  OF  BLOCK  DIMENSION",  15,  /) 

PRINT  6 

6  FORMAT (//IX,  •'  FORWARD  RESIDUALS' ,  60 X,  -  BACKWARD  RESIDUALS  ',  /) 

DO  7  J-l,  NDATM9 

PRINT  8,  (FT  (K,  J),  K-l,  NP),  (BT (K,  U),  K*l,  NP) 

7  CONTINUE 
PRINT  19 

19  FORMAT (//IX,  FORWARD  POWER  BELOW",  60X,  •'  BACKWARD  POWER  BFIOW',/) 

DO  18  I  - 1  >  NP 

PRINT  8,  (P(  I,  K),  K-l,  NP),  (PP(  1,  K.) ,  K=l,  NP) 

18  CONTINUE 

PRINT  3 

3  FORMAT (//IX,  '  FORWARD  REFLECT  COEFF' , 60X,  "  BACKWARD  REFLECT  COEFF 
DO  4  1=1,  NP 

PRINT  8,  (CN( I, K), K=l, NP),  ( CNP ( I . K ) ,  K-“ 1 , NP ) 

4  CONTINUE 

8  FORMAT ( 1 X ,  4HT 2  5,  10X,  *E  12.  5 ) 

PRINT  69 

69  FORMAT (//IX,  '  FORWARD  FILTER  BELOW' , 60X,  'BACKWARD  FILTER  BELOW',/' 

DO  74  K-l,  I9P1 
DO  73  1  =  1 ,  NE* 

PRINT  8,  (CF(  I,  J,  K),  J=l,  NP),  ( CB  ( I ,  J,  K>,  J=l.  NP) 

73  CONTINUE 
WRITE(5,  77) 

71  FORMAT (// ) 

74  CONTINUE 
RETURN 
END 
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SUBROUT I Wb  CM I NV ( A , N,  E ) 

DIMENSION  IPIVOTdM,  A(16,  16),  B(1A>,  INDEX  ( 16,  2) .  PIVOT  (  U) 
DO  20  J=l, N 
20  IF*  I  VOT ( J ) =0 

DO  550  1=1, N 
AMAX=0. 
no  101  J= 1 ,  N 

I F ( I P I VOT  ( J  >  EQ  1 ) GO  T  0  101 
DO  100  K'=  1 ,  N 

I F ( I P I VOT ( K > - 1 > 80 ,  1 00, 7 40 
80  A1 1- ADS ( AMAX ) 

A12=ABS(A( J, K ) > 

IF  (All.  GE.  A 12)  GO  TO  100 
85  I ROW=  J 

ICOLUM=K 
AMAX=A( J,  K) 

100  CONTINUE 

101  CONTINUE 

I P I VOT ( I COLUM >=IPI VOT ( I COLUM ) + 1 
1F<  IRON.  EQ.  I  COLUM)  GO  TO  240 
DO  200  L=1 , N 
SWAP=A<  IRON,  I.) 

A ( I ROW, L ) -A ( I COLUM,  L ) 

200  A ( I COLUM ■ L ) -SWAP 

8WAP=B( I ROW) 

B* ( I ROW ) =6 ( I COLUM ) 

B ( I COLUM  >  =SWAP 
260  INDEX  (1,1  )=-IROW 

INDEX ( I,  2 )  =  1 COLUM 
P I  VOT  <  I )  =A  ( I COI .  UM ,  I  COLUM  > 

A  (  I COIJJM,  I  COLUM )  =  l . 

DO  350  l_=  1 ,  N 

350  A  ( I  COLUM,  L)=A(  I  COLUM,  D/PIVOT  (  I  ) 

B ( I COLUM )=B( I COLUM ) /P I VOT ( I ) 

DO  550  L 1 ~ 1 »  N 
IF  (LI.  EQ.  I  COLUM )  GO  TO  550 
400  T=A(L1, I COLUM) 

A  <1.1,  I  COLUM )  ~0. 

DO  450  L=l,N 

450  A  ( L 1 ,  L )  =A  ( L 1 ,  L )  -  A  ( I  COI  .UM ,  L )  *T 

B ( L 1 ) =B ( L 1 ) -B ( I COLUM ) *  T 
550  CONTINUE 

DO  710  1*1,  N 
L=N+1-I 

IF  ( INDEX  (L,  1).  EQ.  INDEX  (L,  2  >  >00  TO  710 
JROW= I NDF  X ( L  >  1 ) 

JCOLUM- INDEX (L,  2) 

DO  705  K= 1 ,  N 
SWAP*A(K, GROW) 

A ( K ,  JROW ) * A  <  K ,  JCOLUM ) 

A(K, JCOLUM )*SWAP 
705  CONTINUE 

710  CONTINUE 

740  CONTINUE 

RETURN 
END 


SUBROUTINE  STRAND (NOUT.  NR, NDAT, N.  DT,  MRREQ,  NSPE, NPRS) 
DIMENSION  BT(4,  128),ET<4,  12S),CF(4,4,  10),CB(4,  4,  10),FPE(10> 

1.  E<4,  4),  B(4,  4),  0(4,  4),  P(4,  4),  PP(4,  4),  G1  (4,  4),  VI  (4,  4) 

2,  V3  ( 4,  4 ) ,  CN  (4,  4  ) ,  CNP  (4,4),  Ui>  (4,4,  10  ) ,  U6  (  4,  4,  10 ) 

REAL  Q2(4, 4 ) , TI ( 4, 4 > , AA ( 4, 4 ) , BB<4,  4),  CC(4, 4) 

COMMON  INDEX 

COMMON  /BLKO/ET,  BT 
NP 1 - N+ 1 
DO  11  1=1, NP 
DO  11  J=l.  NP 
DO  11  K=1 ,  NP1 
CB(  I,  J.  K )  =0 
CF(  I,  J,  K >=0. 

DO  21  1=1, NP 
CB(  I,  I,  1  )  =  1. 
cf ( r»  i,  i )=i. 

DO  22  1=1, NP 
DO  22  J=l, NP 
P  ( I .  J )  =0. 

DO  21  K= 1 ,  NDAT 

P( I. J)=P( I,  0)+ET (I,  K)*ET(J,  K) 

P ( I , J ) =P ( I »  J ) /FLOAT ( NDAT ) 

PP ( I ,  0 )  =P ( I ,  J ) 

CALL  AKA  IK  (P,  NP,  NDAT.  0.  ,  FF'E,  1 ' 

DO  12  L*2,  N 
LM1-L-1 

NDATM9=NDAT-LM1 
DO  9  K=l, NDATM9 
KP 1 =K+ 1 
DO  9  1=1, NP 
ET ( I ,  K  > =ET ( I , KP 1 ) 

DO  1  1=1, NP 
DO  1  0=1 , NP 
E(I,  J )  =0. 

Bd,  J )  =0. 

0  ( I ,  J )  =0. 

DO  2  K=l, NDATM9 

E<  I,  J)=E<  I,  0)+ET  (I,  K>*ET(  J,  K) 

Bd,  J)  =  =B(I,  J)+BT<  I,  K)*Bf(J,  K> 

0(1,0) -=G (1,0)  +BT ( I,  K) ^E1  < 0,  K ) 

E  ( 1 ,  0 )  =E  (1,0)/  FLOAT  ( NDATM9 ) 

B ( I , 0 ) =B ( I , 0 ) /FLOAT ( NDATM9 ) 

6(1,0)  =-2.  *0  (1,0)  /FLOAT  ( NDATM9 ) 

Q1  ( I»  0)*P(  I,  0) 

Q2  (1,0)  =PP  ( 1 , 0 ) 

CALL  INVERSUH,  NP) 

CALL  I NVFRS ( 02.  NP  > 

CALL  ACCOM  ( NP ,  A  A ,  Q2 ,  B  > 

CAL L  ACCUM ( NP ,  BB , 0 1 ,  F  > 

CALL  ACCUM (NP,  CC,  Q2,  G) 

CALL  GETCNN(  AA,  BH,  CC,  NP,  ON) 

DO  2 3  1=1, NP 
DO  25  0=1, NP 
V1<I,0)=0. 

DO  25  K= 1 , NP 


VI (I,  J)=V1 (I,  J)-CN(K.  I )*PP(K,  J> 

CALL  ACC(  IM  ( NP .  CMP ,  Q 1 .  V 1 ) 

DO  27  1=1, NP 

DO  27  J=1 , NP 

DO  27  K=  j ,  L 

U5<  I.  J,  K)=CF  (  I,  J,  K) 

U6(  I,  J,  K)=CB<  I,  J,  K ) 

DO  4  K=1,L 
INDHX=L— K-l 
DO  4  1=1, NP 
DO  4  NP 

DO  4  KK= 1 , NP 

CF( I, J, K)=CF( I, J, K)+U6< I, KK,  INDEX )*ON(KK,  J> 

CD ( I ,  J,  K)=CE(  I,  J,  K)-*-U5(  I,  KK,  INDEX) ■aCNP { KK ,  J) 

CONTINUE 

DO  5  K=l. NDATM9 

DO  8  1=1,  NP 

V3(I, 1 )=fcT( I,  K> 

V3(  I,  2  )=)?0  (I,  K) 

DO  5  1=1, NP 
DO  5  KK=1,NP 

KT ( I , K ) =HT ( I , K ) +CN ( KK,  I >  *V3 ( KK,  2 ) 

BT ( I , K ) =BT ( I , K ) -CNP <  KK,  I ) *V3 < KK,  1 ) 

CAI  .L  ACCOM  ( NP ,  V 1 ,  PP ,  CN ) 

CALL  ACCOM (NP,  V3,  P,  CNP) 

DO  7  1=1, NP 
DO  7  J=l, NP 
DO  7  K= 1 , NP 

P( I, J)=P( I,  J)-CN(K,  I)*V) (K,  J) 

PP  (I ,  J ) =PP ( I ,  •  J ) -CNP ( K ,  I ) *V3 ( K ,  J ) 

S-FL0A1 (NP*LM1 ) 

CALL  AKA I K ( P , NP , NDAT ,  S , FPE , L ) 

IF  (LM1.  LT.  NPRS ) 00  TO  82 

CAI  L  CPR1NT  (NP,  NDAT ,  I  ..Ml,  P,  PP,  CN,  CNP,  CP,  CD ) 

I F  ( I.M  1 .  LT.  NSPE )  GO  TO  1 2 

CALL  SPEC AL  ( NOOT ,  I  Ml ,  NDAT  ,  NP,  CF,  P,  FPE,  lJ  l  ,  NFRFG!) 

CONTINUE 

RETURN 

END 


SUBROUTINE  GETCNM ( A,  B, C, N, CN) 

REAL  A(4,  4),  B(4,  4),  C(4,  4),  CN(4,  4>,  BB(  1  ) ,  AD (16,  1 6 ) 

N2=N*N 

DO  1  1  =  1,  N 

DO  1  J«1,N 

I J  — ( 1-1 )*N+J 

BB (II) =C ( I ,  J) 

CALL  KRONER  (A,  B,  AB>  N) 

CALL  CMINV(AB,  N2,  BB) 

DO  2  1=1,  N 

DO  2  J=1 ,  N 
Il=( 1-1 )*N+J 
CN( I,  J > =BB (II) 

RETURN 

END 


SUBROUTINE  SPECALtNOUT.  19.  MOAT ,  MP,  OP.  P,  FPL,  DIM  NP HP" Q ) 

REAL  CF(4.  4,  10),  P  <4,  4),  P'PEdO!,  T  1  (  4,  4  ) .  E  ( *■,  4 ) ,  B  ( 4,  4 ) ,  G<  4,  4) 
1,  VI  (4,  4),  VC-  (4.  4),  CU  (4,  4  ),  02(4,  4) 

REAL  SS< 1 6,  123) 

COMMON  ,/BI  Kl/SS 
COMMON  INDEX 
1 9P  1=1 9+1 

ANGLE.=3.  1 4 1 5926/P  LOAT  ( NF  RE  6 ) 

FREQ1  =1.  /  (2.  *  FLOAT  ( NFREQ )  ) 

DEGREE-360.  /(2.  * 3 .  1415926) 

C=l. 

S=0. 

C1=C0S< ANGLE) 

S1=SIN( ANGLE) 

PRINT  2001 

FORMAT  (///I  X.  •"  FREQ  SPEC  11  SPEC22  SPEC33  SPEC44 

1  PH 12  C0H12  PHI 3  C0H13  PH 14  C0H14  PH2S  C0H23 

2  PH24  C0H24  PH34  C0H34  •' , / ) 

DO  99  KK”1,N0UT 

F  REQ=FLOAT ( KK ) *F  REP 1 

C2«=C1*C— S1*S 

S=C1#S+S1#C 

C“C2 

CN  1*1. 

SN1:=0. 

DO  90  I==1.NP 
DO  90  ,J  =1,  IMP 
E  (I,  J>=0. 

B  ( I ,  J )  =0.  , 

DO  91  l.i*l,  I9P1 
DO  92  1=1,  NP 
DO  92  J=1 , NP 

E ( I ,  J )  -:E <  I ,  J)  +CN1*CF <  J,  I ,  L ) 

B ( I ,  J > ~B ( I ,  J )  +SN 1  *CF ( J,  I.L) 

CN2=CN 1 #C-SN 1 *S 

SN1=CN1*S+SN1*C 

CN1=CN2 

CONTINUE 

DO  93  1  =  1,  NP 

DO  93  J=1,NP 

VI  (I,  J)  =  =E<  I,  J) 

CALL  INVP'RS  (VI,  NP ) 

CAI _L  ACCOM  ( NP,  V3,  V 1 ,  B  ) 

CALL  ACCUM(NP, G, B,  V3) 

DO  95  1=1, NP 
DO  95  J=1 , NP 
61(1, J ) =F ( I , J  > +G ( I , J ) 

CALL  INVERS(Q1,  NP) 

CALL  ACCUM(NP, V3, B, VI ) 

CALL  ACCUM ( NP ,  62 ,  Q 1 ,  V3 ) 

DO  96  1  =  1,  NP 
DO  96  J=l, NP 
62 ( I , J ) =-Q2  < I ,  J  > 

CALL  ACCUM ( NP , V 1 , 6 1 , P > 

CALL  ACCUM ( NP, V3, 62,  P) 


DO  30  I  —  1 »  NP' 

DO  30  J~ 1 . NP 
B  ( I .  J )  =0. 

G  ( I  -  J  >  =0. 

DO  30  K=1,NP 

B ( I .  J )  =B <  I .  J ) +V1  ( I .  K  > *Q  1  <  J,  K )  +V3  ( I ,  K )  *02  ( J,  K ) 

30  G( I, J ) =G ( I , J ) +V3 ( I . K)*G1 ( J. K)-V1 (I, K ) *02 ( J , K) 
DO  31  1=1, NP 

DO  31  J=1,NP 
B  ( I ,  J )  =2.  *DT*B  ( I ,  J ) 

31  G  ( I ,  J )  =2.  *DT *6  ( I  *  J  > 

CALL  PHCO(G,  B,  NP,  KK,  FREQ ) 

99  CONTINUE 

CALL  KBP 
RETURN 
END 


SUBROUTINE  PHCO(G,  B, NP, KK, FREQ) 

REAL  G  ( A,  4 ) ,  B  ( 4,  A ) ,  COH ( 6 ) ,  PHAS  ( 6  > ,  SS  (16,  .1  28 ) 

COMMON  /BLK 1 /SS 
DEG=360.  /(2.  *3.  14159265) 

1=0 

NPP=NP-1 
DO  10  J  =1,NPP 
DO  10  K-2, NP 
I F  ( K.  LE.  J )  GO  TO  10 
1  =  1  +  1 

PHAS ( I ) «DEG*ATAN2. ( G  < J,  K ) ,  B ( J, K ) ) 

COH  ( I )  =  ( B  ( J,  K )  **2+G ( J,  K )  **2 )  /  ( B  ( J,  J )  *B  ( K,  K )  ) 

10  CONTINUE 

DO  20  J=1,NP 
20  SS  ( J,  KK )  =B ( J,  J ) 

DO  30  J=l, I 

SS  ( NP+.J,  KK )  "PHAS  ( J ) 

30  SS ( NP+I+J, KK ) =COH ( J ) 

PRINT  2001,  FREQ,  ( B<  1 1 ,  1 1 ) ,  1 1  =  1 ,  NP  > ,  ( PHAS  ( I  > ,  COHO),  1  =  1,  6) 
2001  FORMAT (IX,  F9.  4,  4E9.  2,  12F7.  1  ) 

RETURN 

END _ 

SUBROUTINE  KRONER (A, B, AB, N) 

DIMENSION  A ( 4,  4),  B<4,  4),  AB<  16,  16) 

N2=N*N 
DO  1  1=1, N2 
DO  1  J=1,N2 

1  AB(I,J)=0. 

DO  2  1  =  1,  N 
DO  2  J=1,N 
DO  2  K=1 . N 
11  =  0-1  >*N+K 
I 2= ( J- 1 )*N+K 

2  AB  (Hi  1 2 ) =A ( I ,  J ) 

DO  3  1=1, N 

no  3  J=1 , N 
DO  3  K= 1 , N 
Il=< 1-1 )*N+J 
I2=( 1-1 )*N+K 

3  ABOl,  12 )=AB (II,  I2)+B(K,  J) 

RETURN 

END 


FUNCT I ON  DE  TERM ( ARRAY .  NORDER  > 

D I HENS I ON  ARRAY (4,4) 

DETERM” 1 . 

DO  50  K=l. NORDFR 

IF ( ARRAY ( K.  K).  NE.  0.  >60  TO  41 

DO  23  J  -k  .  NORDER 

I F  ( ARRAY  <  K ,  J ) .  NE.  0.  )  GO  TO  3 1 

CONTINUE 

D£TERI1=0. 

RETURN 

DO  34  I«K. NORDER 
SAVE = ARRAY ( I ,  J ) 

ARRAY ( I . J ) =ARRAY ( I , K ) 

ARRAY ( I , K ) =SAVH 
DETERM--DETERM 
DETERM=»DETERM* ARRAY  ( K,  K) 

IF ( K.  GE.  NORDER) 00  TO  50 
K 1  =  K+ 1 

DO  4 6  I =K1, NORDER 
DO  4<>  J"K1,  NORDER 

ARRAY  ( I ,  J  >  =ARRAY  <  I ,  J )  -ARRAY  ( I ,  K )  *ARRAY  ( K,  .J )  /ARRAY  ( K,  K ) 

CONTINUE 

RETURN 

END _  _  _ 

SUBROUT I NE  I NVERS ( A I NV,  NR ) 

DIMENSION  AINV  ( 4,  4  ) ,  E  INV(  ,  U),BBB<16> 

DO  30  1=1,  NP 
BUB  ( 1  )  =0. 

DO  30  .J=1,NP 
EINVd,  J)=AINV(I,  J) 

C  Al  .L.  CM  I  NV  ( K I N  V ,  NP ,  BOB ) 

DO  31  1*1, NP 
DO  31  J= 1 ,  NP 
AINV(I, J)=EINV( I, J) 

RETURN 

END 


SUBROUTINE  ACCOM <NP, A, B,  C) 
DIMENSION  A(4, 4), B<4, 4 ), C(4, 4) 
DO  15  1*1,  NP 
DO  II  J«1,NP 
A<  I,  J)==0. 

DO  11  K= 1 , NP 

A(I, J)=A< I, J)+B( 1, K)*C(K,  J) 

RETURN 

END 


SUBROUTINE  KBP 
REAL  SS( 5  6,  123) 

COMMON  INDEX 
COMMON  /BLKl/SS 
DO  10  1*1.  16 
INDEX* I 

WRITE< 1 ' INDEX ) (SS( I ,  J) »  J*l,  128) 
CONTINUE 
CAI  L  BELL 
RETURN 


subroutine:  akaikip.  np,  NDAT,  S,  FPE,  INDEX) 
DIMENSION  P<4,  4) ,  FPE< 10), DUM<4, 4) 

DO  1  1=1.  NP 
DO  1  J=1.NP 
DUM ( I •  J)=P<I,  J) 

Dfc TP=  BET ERM ( DUM,  NP) 

FACT  1=FI. OAT  (NDAD+l.  +S 
FACT2=FL0AT ( NDAT ) - 1 .  -S 
F ACT=F ACT 1 /FACT  2 
FPE  ( I NDF  X  )  =DFTP*  ( F ACT**NP  > 

RETURN 

END 
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